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Abstract 

Quantum trajectory methods can be used for a wide range of open quantum 
systems to solve the master equation by unraveling the density operator evolution 
into individual stochastic trajectories in Hilbert space. This C++ class library of- 
fers a choice of integration algorithms for three important unravelings of the master 
equation. Different physical systems are modeled by different Hamiltonians and en- 
vironment operators. The program achieves flexibility and user friendliness, without 
sacrificing execution speed, through the way it represents operators and states in 
Hilbert space. Primary operators, implemented in the form of simple routines acting 
on single degrees of freedom, can be used to build up arbitrarily complex operators 
in product Hilbert spaces with arbitrary numbers of components. Standard alge- 
braic notation is used to build operators and to perform arithmetic operations on 
operators and states. States can be represented in a local moving basis, often lead- 
ing to dramatic savings of computing resources. The state and operator classes are 
very general and can be used independently of the quantum trajectory algorithms. 
Only a rudimentary knowledge of C++ is required to use this package. 
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Program Summary 

Title of program: Quantum trajectory class library 



Program obtainable from: fittp://galisteo.ma.rhbnc.ac.uk/applied/QSD.html| and the au- 
thors. 

Licensing provisions: none 

Operating systems under which the program has been tested: UNIX (Gnu g++), DOS 
(Turbo C++), VMS (DEC C++) 

Programming language used: C++ 

Memory required to execute with typical data: 1MByte 

Has the code been vectorized?: no 

No. of lines in distributed program, including test data, etc.: 8000 

Keywords: open quantum system, master equation, Hilbert space, quantum trajectories, 
unraveling, stochastic simulation, quantum computation, quantum optics, quantum state 
diffusion, quantum jumps, Monte Carlo wavefunction 

Nature of physical problem: 

Open quantum systems, i.e., systems whose interaction with the environment can not 
be neglected, occur in a variety of contexts. Examples are quantum optics, atomic and 
molecular physics, and quantum computers. If the time evolution of the system is ap- 
proximately Markovian, it can be described by a master equation of Lindblad form [1], 
a first order differential equation for the density operator. Solving the master equation 
is the principal purpose of the program. Since the state and operator classes are very 
general, they can be used in any physical problem involving Hilbert spaces with several 
degrees of freedom. 

Method of solution: 

By analogy with the solution of a Fokker-Planck equation by numerical simulation of the 
corresponding stochastic differential equation, a master equation can be solved by simu- 
lating the stochastic evolution of a vector in Hilbert space. The correspondence between 
master equation and stochastic equation is not unique: there are many ways to unravel 
the master equation into stochastic quantum trajectories. The program implements three 
such unravelings, known as the "quantum state diffusion method (QSD)" [2], the "quan- 
tum jump method" [3-5], and the "orthogonal jump method" [6]. The phenomenon of 
phase-space localization [7,8] is exploited numerically by representing quantum states in 
a local moving basis obtained by applying the coherent-state displacement operator to the 
usual harmonic-oscillator basis, often leading to dramatic savings of computing resources. 

Unusual features of the program: 

It is worth emphasizing the effortless way in which operators and states in product Hilbert 
spaces are represented. Primary operators implemented in the form of simple routines 
acting on single degrees of freedom can be used to build up arbitrarily complex operators 
in product Hilbert spaces with arbitrary numbers of components. Building operators, per- 
forming arithmetic operations on operators and states, and applying operators to states is 
done using standard algebraic notation. This program structure has been made possible 
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by systematically implementing object-oriented programming concepts such as inheri- 
tance, concepts which are not (yet) widely used in computational physics. Encapsulation 
of program modules makes it easy to add new basic operators, alternative unravelings of 
the master equation, or different integration algorithms. 

Typical running time: 

The running time depends on the complexity of the problem, the integration time, and 
the number of trajectories required. A typical running time for a simple problem is a few 
minutes. There is no upper limit. 
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Long Write-Up 
1 Introduction 

For many quantum systems of current interest it is no longer possible to neglect the 
interactions with the environment. Those so-called open quantum systems occur in a 
variety of contexts including quantum optics, atomic and molecular physics, and quantum 
computers. Open quantum systems can often be described by a master equation [1|], a 
first-order differential equation for the density operator, in which the internal dynamics 
of the system is represented by the system Hamiltonian H, which is a Hermitian Hilbert- 
space operator, and the interaction with the environment is represented by one or more 
Lindblad operators Lj which are not necessarily Hermitian. 

By analogy with the solution of a Fokker-Planck equation by numerical simulation 
of the corresponding stochastic differential equation (or Langevin equation), a master 
equation can be solved by simulating the stochastic evolution of a vector in Hilbert space. 
The correspondence between master equation and stochastic equation is not unique; there 
are many ways to unravel the master equation into stochastic quantum trajectories [^]. |^. 

The main challenge of this software project was to develop a general program flexible 
enough to accommodate different integration algorithms and unravelings of the master 
equation, as well as the vast range of possible physical systems. In particular, we wanted to 
make it easy to add new algorithms and unravelings, and we wanted a program capable 
of dealing with arbitrary Hamiltonian and Lindblad operators in Hilbert spaces with 
an arbitrary number of degrees of freedom. This task turned out to be ideal for the 
application of object-oriented programming. We chose the C++ language both because 
of its wide availability and because it allowed us to use standard mathematical notation 
for Hilbert-space operations by overloading algebraic operators like '+' and V. 



The core of the program are the C++ classes State and Operator, which represent 
state vectors and operators in Hilbert space. Because of the object-oriented features of 
C++, it is possible to hide the implementation details of these classes completely from the 
classes dealing with the simulation of quantum trajectories. These implementation details 
need not to be known either by a user of the program who wants to choose the quantum 
operators defining the physical problem of interest or by a programmer who wants to add 
a new unraveling of the master equation to the software. A welcome side effect of this 
encapsulation is that the State and Operator classes can be used independently of the 
rest of the code. They should prove useful in many numerical schemes involving Hilbert 
spaces for systems with several degrees of freedom. 

Many Hamiltonian and Lindblad operators can be written as sums of products of 
simple operators acting on a single degree of freedom. Here is an example of a Hamiltonian 
operator coupling a two-level atom (with raising and lowering operators a + and <r_) to 
an electromagnetic field mode (with annihilation and creation operators a and a'): 

H = g (a + a + <7_a^ , (1) 

where the parameter g is the coupling strength. In the following code segment, the atomic 
and field degrees of freedom are labeled and 1, respectively. The Hamiltonian is defined 
in terms of the predefined primary operators SigmaPlus and AnnihilationOperator 
using standard algebraic notation. The class Adapt iveStep is a stepper routine advancing 
the quantum trajectory by a single time step. 

double g = 0.5; 

SigmaPlus Sp(0); // operates on the 1st degree of freedom 

AnnihilationOperator A(l) ; // operates on the 2nd degree of freedom 

Operator Sm = Sp.hcO; // Hermitian conjugate 
Operator Ac = A.hc(); 

Operator H = g*( Sp*A + Sm*Ac ); // Hamiltonian 

AdaptiveStep theStepper ( . . . , H, ...); // ... denotes further arguments 

The important feature illustrated by this example is that the stepper routine is passed 
an object of type Operator without any reference to details like the number of degrees of 
freedom. All the stepper needs to know is that operators can be added, multiplied, etc., 
and that they can be applied to state vectors. 

Internally, the primary operators SigmaPlus and AnnihilationOperator are repre- 
sented as simple loops acting on a single- degree-of-freedom state vector. An instance of 
the more general Operator class is represented by a stack that indicates which primary 
operators are used and the operations by which they are combined. For example, the 
sequence of steps executed by the program when the operator H defined above is applied 
to a state is summarized in the expression 

H\i;)=g(a + (am+^m) , (2) 

in which the elementary steps are applying a primary operator to a state, adding two 
states, and multiplying a state by a scalar. It is clear from this example that a different 
grouping of the terms in the expression for H could lead to inefficient code. This will be 



discussed in Sec. 4.1.2 
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2 Quantum trajectories 
2.1 Master equations 

An open quantum system cannot be described by a Hilbert-space vector evolving 
according to the Schrodinger equation; instead, the state must be described by a density 
operator p whose time evolution generally does not follow any simple law. Fortunately it 
turns out that for a large class of systems the time evolution of the density operator p is 
Markovian to an excellent approximation, i.e., the rate of change of p at time t, dp/dt, 
depends only on p(t), not on the value of p at any earlier time. It has been shown that 
under the Markov approximation the density operator of any open quantum system obeys 
a master equation of Lindblad form |l[ 



where H is the system Hamiltonian and the Lj are the Lindblad operators representing 
the interaction with the environment. 

In many cases, no analytical methods for the solution of the master equation are 
known; one has to use numerical methods. But even a numerical solution of the master 
equation can be very hard. If a state requires D basis vectors in Hilbert space to represent 
it, the corresponding density operator will require D 2 — 1 real numbers; this can often 
be too large for even the most powerful machines to handle, particularly if the system 
involves more than one degree of freedom. 

This problem can be overcome by unraveling the density operator evolution into quan- 
tum trajectories 0, 0, |], ]5[ |], ]7J. Since quantum trajectories represent the system as a 
state vector rather than a density operator, they often have a numerical advantage over 
solving the master equation directly, even though one has to average over many quantum 
trajectories to recover the solution of the master equation. A single quantum trajectory 
can give an excellent, albeit qualitative, picture of a single experimental run. 

2.2 Unravelings 

The three unravelings of the master equation currently implemented are given by the 
following three nonlinear stochastic differential equation for a normalized state vector 




(3) 




l#> = —Hl^dt + ^^ULj-htLj 




(4) 



(ii) the quantum jump equation ||, |^, || 





and (iii) the orthogonal jump equation 0, 0] 



+ £ I - 1 ] 1^}^ . (6) 

The first sum in each of these equations represents the deterministic drift of the state 
vector due to the environment, and the second sum the random fluctuations. Angular 
brackets denote the quantum expectation (G% = {^\G\^) of the operator G in the state 
l^). The d^j are independent complex differential Gaussian random variables satisfying 
the conditions 

Md£j = Md&dtj = , Mdgdtj = Sijdt , (7) 

where M denotes the ensemble mean. The dNj are independent real discrete Poissonian 
random variables satisfying the conditions 

dNf = dN, , dNidNj = , M m dNj = (<LJ4)^ ~ M^jM^iU) dt > ( 8 ) 

where the "conditional mean" M|^\ is defined as the mean over all trajectories for which 
\ip{t)) = \if)), and where A = for the quantum jump equation (||) and A = 1 for the 
orthogonal jump equation @. 

The density operator is given by the mean over the projectors onto the quantum states 
of the ensemble: 

p = M\1>)(1>\- (9) 

If the pure states of the ensemble satisfy one of the quantum trajectory equations (|j), 
(|), or (H), then the density operator satisfies the master equation @: 

Mm))m)\ = p(t), (io) 

where we have assumed that initially the system is in a pure state \ipo) at time t = 0. 
From this it is clear that the expectation value of an operator O is given by 

Tr{6p} = M(ip\6\ip) . (11) 



3 Program Structure 

Our C++ library can be divided roughly into three large parts: 

1. The State class and its associated friend functions. A State includes as member 
data the number of degrees of freedom it represents, how many basis vectors are allocated 
for each degree of freedom, the physical type of each degree of freedom, and (of course) 
the complex amplitudes of each basis vector in the total Hilbert space. The member 
functions include constructors for a number of common State types; arithmetic functions 
enabling States to be added, subtracted, multiplied by scalars, and normalized; functions 
relating to the efficient use of memory, so that a State can be dynamically resized; and 
functions controlling the action of Operators on the State. There are also member 
data and functions relating to the moving basis algorithm, described below. States (and 
Operators) can be used like ordinary variables. In particular, when a locally defined 
State (or Operator) goes out of scope, all memory used by it is properly returned to the 



system; the user of the program need not worry about memory allocation and deallocation 
as this is done automatically. 

2. The Operator class. Operators are defined in terms of their actions on States. 
There is a small class of PrimaryOperators, whose actions on a single degree of free- 
dom are given by pre-defined functions. More complex Operators are defined in terms of 
these PrimaryOperators; they can be added, multiplied, multiplied by scalars or time- 
dependent functions, conjugated, or raised to powers. An Operator's member data in- 
cludes a number of dynamically allocated stacks which indicate which PrimaryOperators 
are used, and the operations by which they are combined. Arithmetic operations on 
Operators are then defined by operations on these stacks. 

3. The Trajectory class and associated classes. These encode the numerical algo- 
rithms for solving the quantum trajectory equations and generating output, with associ- 
ated integration routines, random number generators, and other utilities. Several different 
integration algorithms are currently included, including second- and fourth-order Runge- 
Kutta and Cash-Karp Runge-Kutta with adaptive time steps [§]. These algorithms are 
used to solve the deterministic part of the quantum trajectory equations (^), (^), and (|5[). 
The stochastic terms are solved using first-order Euler integration. The implementation 
of more sophisticated stochastic integration methods (see, e.g., ||) is straightforward. 
Note that it is only in this part of the program that there is any reference at all to the 
details of quantum unravelings. The Operator and State classes are very general. 

These three parts are roughly equal in size, but quite different in internal struc- 
ture. The State class is a single monolithic C++ class with associated functions; the 
Operator class is a parent class with numerous descendent classes representing the differ- 
ent PrimaryOperators. The numerical integration classes are independent of the details 
of State and Operator, and of each other. Because of the object-oriented nature of 
C++, these three groups need know very little about each other's internal workings. 
The following more detailed discussion is not exhaustive; a complete description of the 
code can be found in the extensively commented #include files, particularly in State. h, 
Operator. h, and Traject.h. 

4 The State and Operator Classes 

4.1 One degree of freedom 
4.1.1 States 

We represent a state with a single degree of freedom by an array of N complex 
amplitudes Cj in a given basis {\4>j)}'- 

N 

l^> = (12) 

3=1 

The choice of basis vectors depends on the physical type of the system. For field modes, 
we use Fock states \n); for spins (s = 1/2), we use a z eigenstates | {) and | "[*); for N- 
level atoms, we use energy levels Other types, e.g., molecules or higher spins, can be 
added easily. Of course, a true field mode has an infinite-dimensional Hilbert space. The 
State class represents fields by a finite number of basis states, which should be taken as 
a truncation of the true infinite expansion. 
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To represent a state then requires the physical type (currently FIELD, SPIN or ATOM), 
the number of basis vectors N, and an array of iV complex amplitudes. The state class 
contains constructors for many typical situations. For instance, the expression 

State psi (2, SPIN) ; 

defines psi to be the | j) state of a spin (JV = 2), and 

Complex alpha(0.2,0.3) ; State psi (100, alpha, FIELD) ; 

defines a coherent state |a) with a = 0.2 + 0.3i truncated to N = 100 basis states. 

Arithmetic operations for States are defined internally as operations on the complex 
amplitudes. In the following code examples, the state = 0.5|^i) — \ip2) is formed 
from the Fock states j^i) = |0) and \ip2) = |3), added to \ipi), and then renormalized; 
finally, the inner product z = (V^lV's) i s evaluated. Here N = 10 basis states are more 
than sufficient to represent all states without any truncation. 

State psil(10, 0, FIELD) ; 
State psi2(10,3,FIELD) ; 
State psi3 = 0.5*psil - psi2; 
psil += psi3; 
psi3 .normalize () ; 
Complex z = psi2*psi3; 

The expression psil+=psi3 is superior to the alternative psil=psil+psi3 because it 
avoids the creation of temporary State objects, which is an important consideration in 
high-dimensional Hilbert spaces. 

4.1.2 Operators 

A general way of representing operators is as N x iV complex matrices acting on vectors in 
TV-dimensional Hilbert space. For large N, however, this can be very inefficient, as these 
matrices become very large, and applying them to states requires 0(N 2 ) operations. 
Fortunately, most of the operators of interest in quantum systems are sparse, consisting 
of sums and products of a few primary operators. For FIELDs, such primary operators 
are annihilation and creation operators a and a' and position and momentum operators 
X and P; for SPINs, the primary operators are the Pauli matrices af, for ATOMs, we have 
the transition operators 

In the program, these primary operators are implemented as simple classes, as illus- 
trated for the SPIN operator <j+ in the following code section. 

class SigmaPlus: public PrimaryOperator { 
public : 

SigmaPlus () : PrimaryOperator (0, SPIN) {}; 

SigmaPlus (int freedom) : PrimaryOperator (freedom, SPIN) {}; 

virtual void applyTo (Statefe, int, double) ; 

>; 

void SigmaPlus :: applyTo (Stat e& v, int he, double) { 
switch( he ) { 
case N0_HC: 

v[l] = v[0] ; v[0] =0; break; 
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case HC: 

v[0] = v[l] ; v[l] =0; break; 

} 

> 

The SigmaPlus class is derived from the abstract class PrimaryOperator which serves as 
an interface to the different special classes like SigmaPlus. Apart from the two construc- 
tors, the class contains only the method applyTo. The three arguments of applyTo are 
a single-degree of freedom State, an integer switch determining whether to apply a + or 
its Hermitian conjugate, and a double argument specifying the time for time-dependent 
operators, which is not used here. 

The program represents composite operators, i.e., sums and products of primary oper- 
ators, by stacks containing pointers to primary operators as illustrated in Fig. |l| Those 
stacks are the principal member data of the Operator class, which is the parent class of 
PrimaryOperator and therefore of all special classes derived from PrimaryOperator. For 
a primary operator like SigmaPlus, the stack consists just of the pointer to *this, which 
points to the primary operator itself. Figure |2| shows the hierarchy of operator classes. 

The example stack in box 3 in Fig. [I] is generated by the code segment 

Operator 01 = a + b; 
Operator 02 = (3 * c) * d; 
Operator 03 = 01 - 02; 

where a, b, c, d are assumed to be primary operators defined earlier in the program. 
The example illustrates how addition, subtraction, and multiplication of Operators is 
implemented in terms of operations on the stack. Further operations defined for Operators 
include Hermitian conjugation and raising to an integer power. The C++ inheritance 
mechanism ensures that all these operations are also defined for the derived primary- 
operator classes like SigmaPlus. 

To apply an Operator to a State, the '*' operator can be used as in the following 
example, where psi is a State and 03 is defined above: 

State psil = 03 * psi; 

Internally, this is implemented as a recursive evaluation of the stack. The order in which 
the primary operators are applied in the example can be inferred from the parentheses in 

W> = (a + 6 - 3cd) = (a\tp) + - c (3(d|V») • (13) 

The program keeps the number of operations and the number of temporary States it 
creates to a minimum. Some care has to be exercised, however, to avoid an inefficient 
evaluation order. E.g., in the code segment 

double x=1.5; 
SigmaPlus Sp; 

State psil = 2 . 0*x*Sp*psi ; 

the state Sp*psi is first multiplied by 1.5, then by 2.0, whereas in 

double x=1.5; 
SigmaPlus Sp; 

State psil = (2 . 0*x) *Sp*psi ; 

o 



there is only one multiplication by 3.0. 

The creation of unnecessary temporary States can be avoided by applying Operators 
to States using the '*=' operator as in 

SigmaPlus Sp; 
State psi (2, SPIN) ; 
psi *= Sp; 

When this code segment is executed, no temporary States are created, in contrast to the 
otherwise equivalent code segment 

SigmaPlus Sp; 
State psi (2, SPIN) ; 
psi = Sp*psi; 

A detailed explanation of the stack and the recursive evaluation procedure can be 
found in the extensively commented file Operator, cc. 

4.2 Multiple degrees of freedom 

A quantum system with M degrees of freedom can be represented in a product Hilbert 
space Tii <8> • • • <8> H.m- We assume that there is a finite, perhaps truncated, product basis 
{|0m) <S> • • ' <S> |0n M ) 1 1 — n j — Nj}. Any state \ip) e Hi <S> • • ■ <S> T~Cm can then be written 
in the form 

IV>} = E C ni _ nM \(t> ni )®---®\<t> nM ) , (14) 

ni,...,n M 

requiring a total of N tot = NiN 2 ■ ■ ■ N M basis vectors. To represent a state with multiple 
degrees of freedom, the State class contains as member data the number of freedoms M, 
an array of M physical types, an array of M subspace dimensions Nj, and an array of 
A"tot amplitudes C nij ..., nM . 

Product states can be initialized by passing a list of single- degree-of-freedom states 
to the appropriate State constructor. This is illustrated in the following example, where 
the state |0) <E> |0) <E> 1 1) is assigned to psilni: 

State phi 1(50, FIELD) ; 
State phi2 (50, FIELD) ; 
State phi3(2,SPIN) ; 

State stateList [3] = {phil, phi2, phi3}; 
State psilni(3, stateList) ; 

Entangled states can be constructed by adding several product states or by explicitly 
initializing the array of amplitudes C ni) ... jTlM . 

Operators acting on multiple degrees of freedom are represented as sums and products 
of primary operators each acting on a single degree of freedom. Take the example of a 
primary operator b acting on the first degree of freedom. It can be rewritten as the 
operator b <E> 1 on the product Hilbert space, where 1 is the identity operator acting on 
all the other degrees of freedom. We can write any state as 

IV>> = E IVV...,n M > ® |0n 2 ) ® ' ' • ® \<t>n M ) 5 (15) 
n 2 ,...,n M 

1 n 



the action of 6® 1 on is therefore given by the action of b on the first degree of freedom 
inside a hierarchy of loops over all the other degrees of freedom: 

(b ® 1) |V>> = E (^n 2 ,...,n M )) ® |0n 2 ) ® " • • ® |0n M ) • (16) 
n 2 ,—,TiM 

In the program, the loops are unfolded into one big loop if the primary operator acts on 
the first or last degree of freedom; otherwise the loops are unfolded into two loops, an 
"inner" and an "outer" loop. 

To define, e.g., a primary SigmaPlus operator acting on the n = 3rd degree of freedom, 
the constructor has to be called with the argument n — 1 = 2: 

SigmaPlus Sp(2) ; 

The Operator class is virtually unaffected by the complications arising from multiple 
degrees of freedom (see Fig. |||). Whenever an Operator is applied to a State psi, the 
recursive evaluation of the Operator stack will eventually come across a pointer to some 
primary operator B acting on a particular freedom. At that stage, the pointer to B will 
be passed to the method psi. apply () of the State class, which controls the loops over 
all the other degrees of freedom. Each time the loop is executed, the State class passes 
a single-degree-of-freedom state to the method B.applyToO of the primary operator B. 
The complex amplitudes of this single-freedom state are typically stored at widely spaced 
locations in the array of complex amplitudes C niy __ iTlM , but this fact is completely hidden 
from the primary operator B. 

This way of organizing the program has great advantages. Most importantly, all the 
implementation details of multiple-freedom states are hidden from the Operator class. 
Apart from leading to a transparent program, this makes adding new primary operators 



very easy, as was seen in Sec. |4.1.2| . The definition of the primary-operator class SigmaPlus 
given there is used without modification in the multiple-freedom case. 

Our class library realizes its full potential when all operators are sums and products 
of a few simple primary operators. Although this situation is extremely common in many 
fields, there are important exceptions like the Coulomb potential. While the program 
could be adapted to implement such a case, some of its unique features would be lost in 
the process. 

For efficiency, the State class distinguishes internally between single-freedom and 
multiple-freedom states; many actions are more efficient for a single degree of freedom. 
This distinction, however, is completely transparent. The user need distinguish between 
the two only when constructing the initial state. 



4.3 The Moving Basis 



In quantum-trajectory simulations, one often encounters FIELD states that are well local- 
ized in phase space 
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In cases with strong localization, it is often 
possible to reduce drastically the number N of basis states needed by continually changing 
the basis. If a state is localized about a point (q,p) in phase space far from the origin, 
it requires many number states |n) to represent it. But relatively few displaced number 
states (or excited coherent states) \q,p,n) = D(q,p)\n), are needed, with corresponding 
savings in computer storage space and computation time. The operator D(q,p) is the 
usual coherent state displacement operator JT51, 



D{q,p) 



X P ^ ( pX — qP 



(17) 
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where X and P are the position and momentum operators. The separation of the rep- 
resentation into a classical part (q,p) and a quantum part \q,p,n) is called the moving 
basis [16| or, as in |13|], the mixed representation. To represent a state of type FIELD in 



the moving basis requires to store the complex center of coordinates a — (q + ip)j a/2 in 
addition to the complex amplitudes. A multiple-freedom state in the moving basis with 
several freedoms of type FIELD requires an array of centers of coordinates. 

Implementing the moving basis algorithm is straightforward. Suppose that at time 
t — to the state \ip(t )) is represented in the basis \qo,Po,n), centered at 

(q ,Po) = {{iP{t )m(t )}, (V(to)|P|V>0o)»- (18) 

Then after one discrete time step, the expectations in this basis shift to 

= (Wo + 6t)\X\il>(t + St)), (i/;(t + 6t)\P\il>{t + St))) ^ (q ,p ). (19) 

The computational advantage of a small number of basis states is then retained by chang- 
ing the representation to the shifted basis \qi,Pi,n) centered at q\ and p\. This shift in 
the origin of the basis represents the elementary single step of the moving basis. 

The components of \ij)(tQ + St)) can be computed using the expressions given above. 
The computing time needed for the basis shift is of the same order of magnitude as 
for computing a single discrete time step of one of the quantum trajectory equations. 
Shifting the basis once every discrete time step could therefore double the computing 
time, depending on the complexity of the Hamiltonian and the number of degrees of 
freedom. On the other hand, the reduced number of basis vectors needed to represent 
states in the moving basis can lead to savings far bigger than a factor of 2. 

In the example of second harmonic generation discussed in |T6L two modes of the 



electromagnetic field interact. Using the moving basis reduces the number of basis vectors 
needed by a factor of 100 in each mode. The total number of basis vectors needed is thus 
reduced by a factor of 10000, leading to reduction in computing time by a factor of roughly 
10000/2 = 5000. Furthermore, the fixed basis would exceed the memory capacity of most 
computers. 

The State class includes a variety of basis-changing methods. The most important is 
the method 

void moveCoords ( const Complexfe displacement, int theFreedom, 
double shif tAccuracy ) ; 

which performs a relative shift of the center of coordinates a = (q + ip)/ y/2 by an amount 
given by the complex argument displacement. The integer argument theFreedom spec- 
ifies which degree of freedom is to be shifted — this degree of freedom must be of type 
FIELD. The double argument shif tAccuracy gives the numerical accuracy with which 
to make the shift. The physical state is unchanged by applying moveCoords () , but it 
is represented in a new basis. The method moveCoords () is used in the stochastic inte- 
gration algorithms of the Trajectory class described in Sec. [5| The primary operators 
of type FIELD defined in the files FieldOp.h and FieldOp. cc are implemented in such a 
way that they can handle moving-basis states as well as ordinary states. 

The quantum trajectory equations can contain both localizing and delocalizing terms. 
||. |m |TT[ [l|, 13, |14 . Nonlinear terms in the Hamiltonian tend to spread the wave function 



in phase space, whereas the Lindblad terms often cause it to localize. Accordingly, the 
width of the wave packets varies along a typical trajectory. We use this to reduce the 



computing time even further by dynamically adjusting the number of basis vectors. Our 
criterion for this adjustment depends on parameters e C 1, the cutoff probability, and 
-^pad> the pad size, which represents the number of boundary basis states that are checked 
for significant probability. We require the total probability of the top iV pa d states to be no 
greater than e, increasing or decreasing the number of states actually used accordingly, 
as the integration proceeds along the quantum trajectory. The method of the State class 
used to adjust the basis size is 

void adjustCutof f (int theFreedom, double epsilon, int padSize) ; 

where the arguments specify the degree of freedom to be adjusted, the cutoff probability 
e, and the pad size iVp a d, respectively. 

Like the basis-changing methods discussed above, the method adjustCutof f () is 
typically only used inside integration routines of the Trajectory class. Those methods 
will not normally be called from a top-level program, so the user need not be concerned 
by them. 



5 The Trajectory class 

The Trajectory class and its associated classes, defined in the files Traject.h and 
Traject.cc, implement the integration of the quantum trajectory equations (£|), fl5|), 
and (g). At the heart of this part of the code is the abstract class IntegrationStep 
which serves as an interface for the specific stepper classes implementing single integra- 
tion steps of lenght dt. The stepper classes derived from the class IntegrationStep 
include the class 0rder4Step for a single 4-th order Runge-Kutta step of the QSD equa- 
tion (U) as well as a group of classes using adaptive Cash-Karp Runge-Kutta time steps: 
the class Adapt iveStep for a time step of total length dt of the QSD equation (H), the 
class AdaptiveJump for a time step of total length dt of the quantum jump equation (|5]), 
and the class Adapt iveOrtho Jump for a time step of total length dt of the orthogonal jump 
equation All those classes use a single first order Euler integration step of length dt for 
the stochastic part. Due to the modular structure of the class library, it is straightforward 
to add more sophisticated stochastic integration methods (see, e.g., 0). 



To initialize a stepper, including all temporary memory needed for the integration 
algorithm, all one has to do is call the appropriate constructor as in the code segment 

State psiIni(2,SPIN) ; 
SigmaPlus Sp; 

Operator H = Sp + Sp.hcO; 
int nL = 1 ; 

Operator L [nL] = {0. l*Sp.hc()} 
AdaptiveStep stepper (psilni , H, nL, L) ; 

A less trivial example can be found in the sample program in Sec. Entire quantum 
trajectories are computed by repeatedly calling a stepper from within the Trajectory 
class. A trajectory is initialized as in the following example which is taken from the 
sample program below: 

double dt=0.01; // basic time step passed to the stepper 
ACG gen(38388389) ; // random number generator defined in ACG.h 
ComplexNormal rndm(fegen) ; 



// Gaussian random numbers defined in CmplxRan.h 
Trajectory traj(psilni, dt, stepper, ferndm) ; 

The Trajectory class comprises two methods to launch the simulation, compute ex- 
pectation values of operators of interest, and produce output. The use of the method 
plotExpO, designed to simulate a single trajectory, is explained in Sec. ||. The method 
sumExpO, which is very similar to plotExpO, can be used to compute the mean expec- 
tation values of operators averaged over many trajectories. 



6 Sample program and template 

In this section, we illustrate the main features of the class library in a complete example 
program which can be used as a template. The example program computes expectation 
values for a single trajectory of the quantum state diffusion equation (f|); to compute 
means over many trajectories, one simply replaces the call to traj .plotExpO in the 
template by a call to traj .sumExpO. The system has three degrees of freedom: two 
nonlinear ly coupled field modes described by annihilation operators a\ and a-2, and a 
spin described by raising and lowering operators a + and <r_. The Hamiltonian in the 



interaction picture is [I3J 



H = Ei(a\ — a\) + 7^(ai 2 ci2 — 01^2) + wo+o"- + r]i(a 2 a + — , (20) 

where E is the strength of an external pump field, \ is the strength of the nonlinear 
interaction, u is the detuning between the frequency of the field mode a 2 and the spin 
transition frequency, and rj is the strength of the coupling of the spin to the field mode 
o 2 . The Lindblad operators 

Li = y/^ai , L 2 = ^2r^ 2 a 2 , L 3 = V2kct^ (21) 

describe dissipation of the field modes and the spin with coefficients 71, 72, and k, respec- 
tively. 

The trajectory's initial state is the product state iVw) = |0}(g)|0}(g>| j). The integration 
step-size is dt=0 . 01 and the total integration time is 500*dt = 5. The integration stepper 
AdaptiveStep implements a single time step of length dt of the QSD equation (|j) using 
the Cash-Karp Runge-Kutta algorithm with adaptive time steps ]3J for the deterministic 
part and first-order Euler integration for the stochastic part. 

At times that are integer multiples of 50*dt = 0.5, the expectation values of the 
operators specified in the array out list are computed and written to the files specified in 
the array f list. E.g., the first element of outlist is the operator Xi = <7 + a 2 a_a + . At 
times t = 0, 0.5, . . . , 5.0, the method plotExp computes the expectation values (Ax) and 
var(Ai) = (Ai - (Ai)) and writes t, Re((Ai)), Im((Ai)), Re(var(Ai)), and Im(var(Ai)) 
to the file XI. out. In addition, each time a set of expectation values is computed, the 
program writes 7 numbers to standard output (see the sample output below): the time 
t, 4 expectation values determined by the integer array pipe, the number of basis states 
used, and the number of adaptive steps taken. The integers in the array pipe correspond 
to the columns in the output files containing expectation values (i.e., columns 2 through 
5 of each output file). In the present example, expectation values are computed for the 5 
operators a + a>2cr~a + , a_cr + a2<3"_, a 2 , ni, and h 2 , which are written to 5 output files with 
numbered columns 1 through 20. According to the expression int pipe [] ={ 1 , 5 , 13 , 17}, 
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the expectation values written to standard output are Re((<T + a 2 <J_<j + )), Re((o"_o" + a 2 o"_)), 
Re((ni)), and Re((n2)). 

The moving basis is used for both FIELD degrees of freedom. The basis size is dynam- 
ically adjusted with a cutoff probability e = 0.01 and a pad size N pa d = 2. The sample 
output below shows how the basis size changes with time. Initially, 5000 = 50 * 50 * 2 
states are allocated, but at time t = 0.5, only 18 states are needed. Subsequently, the 
basis size fluctuates around a typical size of 70 states. 

Here is the complete program: 

#include "Complex. h" 

#include "ACG.h" 

#include "CmplxRan.h" 

#include "State. h" 

#include " Operator. h" 

#include "FieldOp.h" 

#include "SpinOp.h" 

#include "Traj ect . h" 

int main() 
{ 

// Primary Operators 

AnnihilationOperator A1(0); // 1st freedom 
NumberOperator N1(0); 

AnnihilationOperator A2(l); // 2nd freedom 
NumberOperator N2(l); 

SigmaPlus Sp(2); // 3rd freedom 

Operator Sm = Sp.hcO; // Hermitian conjugate 

Operator Acl = Al.hcO; 

Operator Ac2 = A2.hc(); 
// Hamiltonian 

double E = 20.0; 

double chi = 0.4; 

double omega = -0.7; 

double eta = 0.001; 

Complex 1(0.0,1.0) ; 

Operator H = (E*I) * (Acl-Al) 

+ (0.5*chi*I)*(Acl*Acl*A2 - Al*Al*Ac2) 
+ omega*Sp*Sm + (eta*I) * (A2*Sp-Ac2*Sm) ; 
// Lindblad operators 

double gammal = 1.0; 

double gamma2 = 1.0; 

double kappa = 0.1; 

const int nL = 3; 

Operator L [nL] ={sqrt (2*gammal) *A1 , sqrt (2*gamma2) *A2 , sqrt (2*kappa) *Sm} ; 
// Initial state 

State phil (50, FIELD) ; // see Section 4.2 

State phi2 (50, FIELD) ; 
State phi3(2,SPIN) ; 

State stateList[3] = {phil,phi2,phi3}; 

1 c; 



State psiIni(3,stateList) ; 
// Trajectory 

double dt = 0.01; // basic time step 

int numdts =50; // time interval between outputs = numdts*dt 

int numsteps =10; // total integration time = numsteps*numdts*dt 
int nOf MovingFreedoms = 2; 

double epsilon =0.01; // cutoff probability 

int nPad =2; // pad size 

ACG gen(38388389) ; // random number generator with seed 

ComplexNormal rndm(fegen) ; // Complex Gaussian random numbers 
AdaptiveStep stepper (psilni , H, nL, L) ; // see Section 5 

Trajectory traj (psilni, dt, stepper, ferndm) ; // see Section 5 
// Output 

const int nOfOut = 5; 

Operator out list [nOf Out] ={ Sp*A2*Sm*Sp, Sm*Sp*A2*Sm, A2, Nl, N2 }; 
char *f list [nOf Out] ={"X1 . out" , "X2 . out" , "A2 . out" , "Nl . out" , "N2 . out "} ; 
int pipe [] = { 1, 5, 13, 17 }; // controls standard output 
// Simulate one trajectory 

traj .plotExp( nOfOut, outlist, flist, pipe, numdts, numsteps, 
nOf MovingFreedoms , epsilon, nPad ); 

} 

In addition to the output files XI . out, X2 . out, A2 . out, Nl . out, and N2 . out, the program 
writes the following lines to standard output: 




























5000 
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0, 


,000505736 


0. 


.000504849 


52. 


.3875 


3, 


.5807 
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1 




0, 


.0131402 


0. 
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.8747 
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.1089 
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1. 


,5 


0. 


.0329714 


0. 


,0320222 


32. 


.8707 


44. 


.3184 
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50 


2 




0. 


. 0425276 


0. 


.0455457 


32. 


.1562 


41. 


.7798 


70 


56 


2. 


.5 


0. 


.0284912 


0. 


,0564117 


34. 


.85 


37. 


.8809 


80 


117 


3 




0, 


.0260639 


0. 


.0626976 


33. 


.9828 


39. 


.3437 


80 
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3. 


.5 


0. 


. 0544306 


0. 


,0439029 


51. 


.0632 


37. 


.6462 


70 


99 
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0, 


.0796275 


-0. 


.0209383 


41. 


.9614 


38, 


.0884 


70 


167 


4, 


.5 


0, 


.0834672 


-0. 


,0543796 


33. 


.1194 


36. 


.1007 


70 
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5 




-0, 


.00616844 


0. 


.0110794 


76. 


.6321 


29. 


.4303 


50 


119 
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Figure 1: These examples show how the internal stack representations of primary and 
composite Operators are combined in arithmetic operations. Notice that while arithmetic 
expressions are parsed from left to right, the order in which Operators are applied to 
States is from right to left. 



Figure 2: In this diagram, the arrows point from parent classes to derived classes. The 
classes listed in each box are declared in the #include file given above the box. Arithmetic 
operations are defined in the Operator class. The PrimaryOperator class serves as an 
interface for the specific FIELD, SPIN, and ATOM operators. Adding operators of either an 
existing or a new type is straightforward. 



Figure 3: When a pointer to a primary operator acting on a particular freedom is encoun- 
tered during the evaluation of an Operator stack, control is passed to the State class, 
where within loops over the basis states of all the other degrees of freedom, the primary 
operator is applied to a succession of single-freedom states. This means that the Operator 
class and its derived classes do not need to distinguish between single and multiple free- 
dom states; all details concerning the multiple-freedom case are hidden within the State 
class. 
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First declare the primary operators 
a, b, c, d; 



1 

2 
3 



Operator 01 
Operator 02 
Operator 03 



= a + b; 

= (3 * c) * d; 

= 01-02; 



Operator.h 



Operator 

+ — * += 

hc(), powC 



PrimOp.h 



FieldOp.h 



AnnihilationOperator 

a = (x + ip)/y/2 
NumberOperator 

n = a^a 
IdentityOperator 

1 

XOperator 

x = (a + a x )l\l2 
POperator 

p = i(a^ — a) I y/2 



PrimaryOperator 

applyTo () 



SpinOp.h 



SigmaX 



(Tr, 



SigmaY 



SigmaZ 



SigmaPlus 



AtomOp.h 




Operator 



Stack points 
to primary 
operators which 
act on particular 
freedoms. 



Apply each primary operator to the state; 
tranformed state returned. 



apply () 



applyTo () 



PrimaryOperator 



Loop over all other degrees of freedom. 



State 

M degrees of freedom. 



Pass succession of one-freedom states 
to primary operator; transformed states 
returned. 



Act on single freedom states. 



